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The transition from a liquid to a glass remains one of the most poorly understood phenomena in condensed 
matter physics, and still no fully microscopic theory exists that can describe the dynamics of supercooled liquids 
in a quantitative manner over all relevant time scales. Here we present a theoretical framework that yields 
near-quantitative accuracy for the time-dependent correlation functions of a glass-forming system over a broad 
density range. Our approach requires only simple static structural information as input and is based entirely 
on first principles. Owing to its ab initio nature, the framework offers a unique platform to study the relation 
between structure and dynamics in glass-forming matter, and paves the way towards a systematically correctable 
and ultimately fully quantitative theory of microscopic glassy dynamics. 


Understanding the dynamics of supercooled liquids repre¬ 
sents one of the major challenges in condensed matter science 
(2]. Perhaps the most striking feature of vitrification is the 
observed dramatic increase in viscosity (or relaxation time) 
upon only a relatively mild change in thermodynamic control 
parameters, e.g. temperature or density fi. This highly non¬ 
linear response of the system is further accompanied by only 
subtle changes in the microscopic structure, thus raising the 
question as to what physical mechanism underlies the glass 
transition fill. 


Among the various theories of the glass transition proposed 
in the last few decades fi], mode-coupling theory (MCT) has 
acquired a prominent place in this field of research |6 Kj. 
MCT is the only strictly first-principles theory of glassy dy¬ 
namics, and it can accurately predict many features of the 
time-dependent dynamics of supercooled liquids from sim¬ 
ple static information alone, including multi-step relaxation 
patterns fi, stretched exponentials fi], and growing dynami¬ 
cal length scales fi]. MCT forms a mean-field framework for 
broader theories such as Random First-Order Theory (RFOT) 
0, but it remains a major open question how the relaxation 
scenarios posited in RFOT can be captured from first princi¬ 
ples and fully microscopically along the same lines as MCT. 
Thus far, efforts to systematically correct the MCT e quat ions, 
e.g. via higher order field-theoretic loop expansions II lUl41. 
have not been explicitly carried out in a practical setting. 
Hence, unlike the situation that arises in fields such as elec¬ 
tronic structure theory, where the correlation energy can be 
systematically captured by building upon Hartree-Fock the¬ 
ory m, there exists no means of providing a more accurate 
description of glassy dynamics by microscopically correcting 
the lowest-level mean-field approach. Indeed, a commonly 
espoused viewpoint is that there exists a ”no-theory region" 
between MCT and the true glassy regime 0. In this work 
we show that such a theory can in fact be formulated and suc- 
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cessfully carried out, and we present comparisons with dy¬ 
namics in a realistic model of a glassy hard-sphere suspension 
to quantify the systematic improvement achieved with our ap¬ 
proach. 

The key quantity in our discussion is the two-point density 
correlation function F(k,t) = N~ 1 (p^(0)py i (t)), which 
probes correlated particle motion over time t and at inverse 
length scale k. Here Pk(t) represents a density mode at time 
t and wavevector k , k = |k|, TV is the total number of par¬ 
ticles, and the brackets denote a canonical ensemble average. 
The exact time evolution of F(k, t) is governed by an integro- 
differential equation with memory function M(k,t), the lat¬ 
ter of which accounts for time-dependent damping effects fi. 
More precisely, M t) is related to the autocorrelation func¬ 
tion of a fluctuating force, which captures the effect of all de¬ 
grees of freedom orthogonal to the density field fi. It can 
be shown that this fluctuating force is dominated by pairs of 
density modes PqPk-q, allowing the memory function to be 
expressed in terms of/rwr-point density correlation functions 
(p_q(0)p q _k(0)pq'(t)pk-q'(i)) fi Q. Within the tradi¬ 
tional MCT framework, the four-point correlation functions 
are subsequently approximated as products of two two-point 
correlation functions, thereby rendering a closed, non-linear 
equation of motion for F(k , t) fifi. The microscopic struc¬ 
ture of the system enters through the static structure factor 
S(k) = F{k, 0) 0, after Gaussian and convolution approx¬ 
imations are made for the statics. 

In the present work, we seek to avoid the uncontrolled 
‘Gaussian’ factorization of the dynamic four-point correla¬ 
tions, instead retaining these correlations in the expression 
for F(k,t) and developing a new exact equation of motion 
for the memory term. As shown explicitly in the Supple¬ 
mentary Information (SI) ififi . the dynamics of these four- 
point correlation functions is governed by a new memory 
kernel that (after projection onto triplet-density modes) con¬ 
tains six-point correlations, i.e. correlation functions of the 
form (pkAOjpkzWpks^Pk^Pk^p^t)). These are 
controlled by eight-point correlations, and so on, thus allow¬ 
ing one to delay the uncontrolled factorization approximation 
to a later stage. Our approach builds upon the important work 
of Szamel, who first introduced such a hierarchical scheme 
of microscopic kinetic equations 1811 . Several recent stud- 
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ies have demonstrated that this approach indeed holds great 
potential as a microscopic theory of glassy dynamics. In par¬ 
ticular, it was shown that higher-order correlations can sys¬ 


tematically improve the predicted MCT critical point 181 1 


Moreover, two schematic (fc-independent) proof-of-principle 
studies predicted a wealth of novel relaxation patterns for con¬ 
tinuous, discontinuous, and avoided glass transitions MM- 
For example, the theory may describe both Arrhenius and 
super-Arrhenius growth of the relaxation time 112111 , thus pro¬ 
viding potentially the first rigorous first-principles rationale 
for the concept of fragility |2, 22h . 


Here we report on the first fully microscopic, high-order 
calculations to describe the explicit time- and /,'-dependent dy¬ 
namics of a realistic glassy system. We focus on a system 
composed of quasi-hard spheres, and compare our theoreti¬ 
cal predictions with the exact dynamics obtained from com¬ 
puter simulations ll23il . By using only the simplest measure 
of static correlations, i.e. S(k ), as sole input to our theory, we 
already achieve high quantitative accuracy for F(k. t), for all 
wavevectors, in the regime from low to moderate supercool¬ 
ing. In view of its first-principles nature, the theory can shed 
important new light on the fundamental connection between 
structure and dynamics in glass-forming matter, and paves 
the way towards an ultimately fully quantitative description 
of glassy dynamics. 


The first main result of this work is the formulation of a 
new set of microscopic equations of motion for the dynamic 
density correlation functions of a glass-forming system. The 
full derivation is described in the SI. We consider the dynam¬ 
ics of the normalized 2n-point density correlation functions 
... ,k n ,t), which probe particle correlations over n 
distinct fc-values. 
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Within our framework, these correlators obey the general 
equations of motion 


,...,k n ,t) + u^ n (ki ,..., k n , t) 

,..., k n )$ n (k 1 ,...,k n ,t) 

+ [ M n (ki,... ,k n ,T)$ n (ki,... ,k n ,t - r)dr = 0, (2) 
Jo 

where the dots denote time derivatives, v represents a coeffi¬ 
cient that accounts for the short-time dynamics, and the bare 
frequencies are given by 
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Here k B is the Boltzmann constant, T is the temperature, and 


m is the particle mass. For the memory functions we have 

pk B T ^ nl(ki) 
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where p is the total density, S,_ :) is the Kronecker delta func¬ 
tion, and fqki-q are static vertices that represent wavevector- 
dependent coupling strengths for the higher-order correla¬ 
tions. These vertices depend on the microscopic structure of 
the system contained in S(k), which serves as input to the 
theory fH. Note that the memory function at level n con¬ 
tains n different terms, each of which contains a 2(n + 1)- 
point correlator that measures correlations over a distinct set 
of wavenumbers {q, Aq,..., |ki — q|, ki + i, ..., k n }. In ar¬ 
riving at the above equations, we have retained all diagonal 
dynamic 2n-point correlation functions, and employed Gaus¬ 
sian and convolution approximations for the static correlators 
(see SI). These remaining (uncontrolled) assumptions are es¬ 
sentially the same as at the standard-MCT level; the crucial 
difference between our work and standard MCT is that we 
do not factorize the diagonal dynamic multi-point correlators. 
We expect the inclusion of off-diagonal dynamic correlators 
to lead to a matrix form of the hierarchy that is likely to be 
far more demanding computationally, but identical in over¬ 
all structure compared to the present set of equations. We 
may close this hierarchy of coupled equations at an arbitrary 
level n = N by applying a ‘mean-field’ (MF) factorization 
of the form ~ [the simplest example 

of which is the standard-MCT closure $ 2 (f) ~ 4 > i(f)<I > i(f)] 
or by simply setting &N{t) = 0. The latter implies that 
~ exp(— VL 2 N _ l t), hence we refer to this trunca¬ 
tion as an exponential closure. Finally, we emphasize that 
our theory is free of any fitting or rescaling parameters, thus 
providing dynamic predictions from purely static information 
alone. 

In Fig. [T] we show our numerical results for the density 
correlator $1 (k,t) = F(k , t)/S(k) for the quasi-hard sphere 
system of Ref. [23]. While this system is slightly polydis- 
perse to prevent crystallization, we have used the averaged 
(effectively monodisperse) structural properties as input to our 
calculations, as was also employed by the authors of Ref. 

1230. As a stringent test of our theory, we compare the the¬ 
oretically predicted dynamics with the numerically simulated 
<1» | (k. t ) for different wavevectors and different volume frac¬ 
tions (f> = pnd 3 /6, where d is the average particle diame¬ 
ter. Note that, while our theory also readily gives access to 
the higher-order correlators, the simulated data are limited 
to The wavevectors we consider correspond to the 

main (kd = 7.4) and second (kd « 13.0) peak of the static 
structure factor. We have solved the time-dependent coupled 
equations numerically employing both MF and exponential 
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FIG. 1: Two-point density correlation functions <f>i (k,t) for quasi- 
hard spheres at different volume fractions and wavevectors: (a), <j> = 
0.450, kd = 7.4, (b), <f> = 0.500, kd = 7.4, (c), </> = 0.570, 
kd = 7.4, (d), <f) = 0.570, kd = 13.0. The theoretical results have 
been obtained using either exponential (red lines) or MF (blue lines) 
closures, with closure levels N = 2,3, and 4. The N = 2 MF 
closure is equivalent to standard MCT. All legends are as in panel 
(a). The simulated data (black circles) are taken from Ref. |23i]. 


closures for truncation levels N = 2,3,4. Computational de¬ 
tails can be found in the SI. 

Before discussing Fig. [T] in detail, we first point out an im¬ 
portant general observation with regard to the convergence of 
our results. For all volume fractions and wavevectors consid¬ 
ered, we find that the MF closures provide an upper bound to 
the dynamics, i.e. the N = 2 (standard-MCT) closure pre¬ 
dicts the slowest relaxation, while the exponential closures 
give a lower bound to the dynamics. This is expected because 
any mean-field theory will generally underestimate the ef¬ 
fects of ergodicity-restoring fluctuations, consequently over¬ 
estimating the dynamic slowing down. Conversely, an ex¬ 
ponential closure rigorously ignores higher-order memory ef¬ 
fects that may drive the system into the glassy state [since 
~ Mjy + i(f) = 0 ], resulting in relaxation patterns that 
are too fast. With increasing closure level N, however, the two 
types of closures systematically converge towards each other 
over an increasingly large time domain. This uniform con¬ 
vergence pattern can be rigorously demonstrated within the 
confines of a schematic model |[ 20 ], and provides a unique 
and clearly achievable notion of convergence even though no 
small parameter exists in our approach. 

Let us now focus on the time-dependent relaxation behavior 
of the quasi-hard sphere system for specific values of </>. For 
the lowest volume fraction considered, <f> = 0.450, the system 
is still strongly ergodic and 'tq (/,:, t) exhibits relatively simple 
relaxation [Fig. [TJa)] . Here our N = 4 data are almost per¬ 
fectly converged upon the simulated 4>i (k,t) curve over the 
complete time domain. Note that the standard-MCT result is 
also reasonably accurate, but slightly overestimates the molec¬ 


ular relaxation time. These findings hold for all wavevectors 
considered. Thus, in the ’’normal liquid" regime, the predic¬ 
tions of our higher-order theory are virtually exact with re¬ 
spect to computer simulations. 

For a higher volume fraction, <j> = 0.500, the relaxation 
time is seen to increase and the system enters the supercooled 
regime [Fig. [ljb)]. Again we find nearly perfect quantita¬ 
tive agreement between our N = 4 results and the simulated 
density correlator at kd = 7.4, indicating that our higher- 
order microscopic approach accurately captures the onset of 
glassy dynamics at the correct absolute value of </>. This is 
to be contrasted with standard MCT (N = 2), which clearly 
overestimates the glassiness of the system and can only re¬ 
produce the simulated results for a rescaled volume fraction 
<f >mct < 0.500 J 23 I . Since the higher-order theory presented 
here does not require such rescaling of the control parame¬ 
ter, the inclusion of higher-order dynamic correlations thus 
appears to capture all relevant relaxation mechanisms in the 
weakly supercooled regime. 

We now turn our attention to the highest volume fraction, 
f = 0.570, where the system is significantly supercooled. 
Standard MCT predicts a diverging relaxation time above 
^mct _ 0 500 f 23 ll . as seen j n our N — 2 results in Figs. 
[Tfc)-(d). The simulated <h 1 (k. t ) curves, however, show com¬ 
plete relaxation to zero, indicating that the system is still er¬ 
godic at this density. In agreement with these simulations, our 
N = 3 and N = 4 results also predict full relaxation, thus 
restoring ergodicity and rounding off the spurious MCT tran¬ 
sition. We find good agreement with the simulated data for 
both wavevectors, but our higher-order theoretical results ap¬ 
pear to converge toward a relaxation time that is slightly too 
fast compared to the simulations. The predicted curves for 
kd = 13.0 also develop a small shoulder at intermediate times 
that is absent in the simulated data. This slightly poorer agree- 
ment for kd > 7.4 was also noted by Weysser et al. 112311 for 
the standard-MCT case, even after rescaling of the data. The 
origin of the discrepancy may be related to the fact that the 
use of a density-mode basis neglects some short-time correla¬ 
tions, which are expected to be more pronounced away from 
the first peak of S(k). We attribute further discrepancies to 
the remaining approximations in our theory, i.e. the neglect of 
off-diagonal dynamic correlation functions (p q (0 )... p q /(f)) 
with q ^ q', and the factorization of the static multi-point cor¬ 
relations. Indeed, one intriguing possibility is that the weaker 
agreement with simulation at the highest densities indicates 
that multi-point static correlations, such as the much studied 
point-to-set correlations J 3 I], begin to extend a substantial in¬ 
fluence on the relaxation. In this regard, our approach also 
has utility in delineating the role of static structural properties 
beyond S(k) on the dynamics. Finally, Weysser et al. have 
found that the standard-MCT predictions for this system can 
be improved, in particular for small /.--values, by accounting 
for polydispersity effects 112311 . Our theoretical results may 
thus also be further improved by careful treatment of the sys¬ 
tem’s polydispersity. 

Despite the quantitative differences found in Fig. |TJc)— (d). 
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FIG. 2: Root-mean-square deviations between the four-point density 
correlator $2 (fci, fe, t ), as obtained from our high-order calculations 
under an TV = 4 MF closure, and the standard-MCT approximation 
4>^ CT (fci, t)4>5 /ICT (fc2, t), as a function of wavevectors ki and 
at (a) < j> = 0.450 and (b) f = 0.570. 


it is evident that our higher-order framework in its current 
form captures a significant portion of the activated dynamics 
beyond the standard-MCT scenario, at least up to moderate 
volume fractions. We emphasize that our theoretical predic¬ 
tions are free from any fitting or rescaling procedures, and that 
regular MCT would not be able to achieve a similar degree of 
accuracy, even after introducing a rescaling parameter. The 
good agreement we find between our work and simulations is 
highly non-trivial, considering that the only input to our the¬ 
ory is S(k), the simplest measure of the microscopic structure 
of the system. Our results thus support the view that acti¬ 
vated dynamics-at least in the weak to moderate regime of 
supercooling-can be described via a strictly dynamic frame¬ 
work, in the form of deeper memory effects that occur over 
increasingly many wavelengths. 

As a final part of our analysis, let us elaborate on how our 
theory improves upon the standard-MCT framework in terms 
of the memory function. Recall that standard MCT approxi¬ 
mates the lowest-order memory contribution <I> 2 (<7, |k — q|. t) 
as ^ CT (g,|k-q|,t) ~ $r CT (g,f)^ CT (|k-q|,i) 

|8(]. Conversely, our framework retains the four-point 
correlations explicitly and treats their dynamics rigorously 
through a hierarchy of equations. Understanding how the 
<T» 2 (q, |k — q|, t) that emerges from our approach differs from 
qjMC'r jjj. _ q|,t) thus provides insight into the improve¬ 
ments afforded over standard MCT. In Fig. [2] we com¬ 
pare our best estimate ( N = 4) for the dynamics of 
the four-point correlation function fa, t) with the 

standard-MCT approximation $^ CT (fci, f)$5 4CT (fc2, f). We 
plot the root-mean-square deviation (RMSD) |$ 2 (fci, k 2 ,t) — 
jjjMCT ^ q>MCT , f) | / \/Nt, where Nt is the total num¬ 
ber of time points, for volume fractions <f> = 0.450 and 
<fi = 0.570. It may be seen that for both volume fractions 
the RMSD is modulated by the structure of S(k). In the low- 
density regime = 0.450, Fig. [2ja)], the largest deviations 
are found at large wavevectors. These differences arise from 
the interplay between the bare frequency term D 2 [k\ , kf) and 
the inertial term $2 (k \, k 2 , t ), which make the four-point cor¬ 
relation functions $ 2 (k±,kz,t) oscillate and become negative 
at short times. Such negative contributions are mostly absent 


in standard MCT, since the diagonal MCT terms k, t) ~ 
|4>i(fc,f )| 2 are always positive. Hence, standard MCT tends 
to overestimate the memory term M) (k. t) and consequently 
overestimates the relaxation time. These oscillatory patterns 
only arise in the underdamped limit, however, and neglecting 
the inertial term (i.e. overdamped dynamics) dramatically re¬ 
duces the RMSD in the large-wavevector regime. 

In the deeply supercooled region = 0.570, Fig. HJb)], 
the overall RMSD increases, and the largest deviation is now 
observed at smaller wavevectors, most notably near the main 
peak of S(k) at kid « k 2 d « 7.4. The RMSD also exhibits 
local maxima at wavevectors of kd ~ 13.0, corresponding to 
the second peak of the structure factor, and kd ss 1.0. Thus, 
the ergodicity-restoring processes that round off the standard- 
MCT transition are mainly contained in the lower-/,: behav¬ 
ior of the memory function. It is well established that the 
standard-MCT framework is generally the least accurate in 
the hydrodynamic k —> 0 limit, resulting in e.g. an underesti¬ 
mation of Stokes-Einstein violation in the deeply supercooled 
regime 1241] . Since the higher-order corrections at large vol¬ 
ume fractions occur mostly at small fc-values, one may ex¬ 
pect that our theory might improve the predicted breakdown of 
the Stokes-Einstein relation and other related quantities. This 
possibility will be tested in future work. 

In conclusion, we have presented a novel first-principles 
theory to describe the microscopic dynamics of glassy ma¬ 
terials with near-quantitative accuracy in the low to moder¬ 
ately supercooled regime. The framework requires only static 
structure as input and has no free parameters. Our results 
demonstrate that a quantitative account of dynamical corre¬ 
lations may be generated from only two-body statics-even in 
a regime where many-body statics and dynamical heterogene¬ 
ity are markedly growing-, thus shedding light on the connec¬ 
tion between structure and dynamics in supercooled liquids. 
Our framework may also eventually be used to provide pre¬ 
dictive insights into regimes of glassy behavior inaccessible to 
molecular-dynamics simulations. Not only does our method 
give straightforward access to e.g. long times and higher-order 
correlation functions, it also suggests that by incorporating in¬ 
creasingly many (static and dynamic) correlations, one may 
achieve systematically improvable and ultimately fully quan¬ 
titative accuracy in the deeply supercooled regime. 

As a final note, it should be remarked that the hierar¬ 
chy presented here resembles the Martin-Schwinger hierarchy 
1 25112611 in field theory, where low-order Green’s functions are 
hierarchically connected to higher-order ones through an infi¬ 
nite set of coupled equations. In the context of the present 
work, we point out that there has been substantial recent 
progress in the formulation of field-theoretic approaches for 
the glass transition that allow for a systematic inclusion of 
higher-order terms that are omitted from standard-MCT-like 
theories |27|-29[]. It is likely that such approaches can be for¬ 
mulated in a Martin-Schwinger-type form, which would en¬ 
able the direct comparison of the hierarchy presented here 
with recent field-theoretic formulations. Such a comparison 
might also make clearer the effect of several approximations 
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that have been made here. Future work should be devoted to 
such a program. 
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